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Abstract:  The  effect  of  various  impurities  and  micro-alloying  additions  (B,  N,  C, 
O,  Al,  Si,  S  and  P)  on  the  intrinsic  resistance  of  the  23  (111)  grain  boundary  in 
tungsten  has  been  investigated  using  the  molecular  dynamics  simulation.  The 
atomic  interactions  have  been  accounted  for  through  the  use  of  Finnis-Sinclair 
interatomic  potentials.  The  fracture  resistance  of  the  grain  boundary  has  been 
characterized  by  computing,  in  each  case,  the  ideal  work  of  grain  boundary 
separation,  the  mode  I  stress  intensity  factor  and  the  Eshelby’s  F,  conservation 
integral  at  the  onset  of  crack  propagation.  The  results  obtained  suggest  that  pure 
tungsten  is  relatively  resistant  to  grain  boundary  decohesion  and  that  this 
resistance  is  further  enhanced  by  the  presence  of  B,  C  and  N.  Elements  such  as 
O,  Al  and  Si  however,  have  a  relatively  minor  effect  on  the  cohesion  strength  of 
the  S3  (111)  grain  boundary.  In  sharp  contrast,  S  and  P  greatly  reduce  this 
strength  making  tungsten  quite  brittle.  These  findings  have  been  correlated  with 
the  effect  of  the  impurity  atoms  on  material  evolution  at  the  crack  tip.  ©  1997 
Elsevier  Science  Ltd.  All  rights  reserved 


INTRODUCTION  density  of  THAs  is  only  slightly  lower  than  that 

of  pure  tungsten,  the  alloys’  anti-armor  per- 
Tungsten-base  alloys,  due  to  their  high  density  formance  is  quite  inferior  relative  to  that  of 

and  strength,  are  widely  used  in  army  systems,  depleted  uranium,  a  competing  anti-armor 

especially  for  anti-armor  applications.  However,  material.  This  difference  is  related  to  the  fact 

technically  pure  tungsten  is  extremely  brittle  that,  while  a  depleted  uranium  penetrator 

and  its  ductile-* brittle  transition  temperature  develops  adiabatic  shear  bands  upon  impact 

(DBTT)  is  typically  as  high  as  300°C.  To  over-  which  give  rise  to  a  ‘self-sharpening’  effect2,  the 

come  the  problems  associated  with  the  limited  THA  penetrators  ‘mushroom’,  which  leads  to 

toughness  of  tungsten,  the  so-called  tungsten  an  increase  in  the  penetrator  diameter  and,  in 

heavy  alloys  (THAs)  have  been  developed1.  turn,  to  a  reduction  in  the  depth  of  penetration. 

Manufactured  by  liquid  phase  sintering,  these  Owing  to  the  well-known  environmental  prob- 

alloys  are  actually  metal-matrix  composites  lems  associated  with  using  depleted  uranium, 

consisting  of  hard  tungsten  particles  embedded  significant  resources  are  currently  being 

in  a  relatively  mild  Ni-Fe  matrix.  Although  the  directed  towards  developing  THAs  with 
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improved  adiabatic  shear  band  behavior.  The 
usual  approach  is  to  significantly  modify,  or  find 
an  alternative  to,  the  Ni-Fe  matrix. 

A  different  approach  has  recently  been  initi¬ 
ated  by  Krasko3,  who  suggested  micro-alloying 
as  a  way  of  obtaining  a  ‘ductile’  tungsten.  Being 
a  B.C.C.  metal,  such  tungsten  is  expected  to 
acquire  the  necessary  adiabatic  shear  band 
behavior  since  the  latter  has  already  been 
observed  in  high-strength  B.C.C.  steels  in  condi¬ 
tions  of  ‘marginal’  ductility4. 

The  reduced  cohesion  of  grain  boundaries  is 
frequently  cited  as  a  major  factor  limiting  duc¬ 
tility  and,  in  turn,  the  performance  and  relia¬ 
bility  of  high-strength  metallic  materials5,6. 
Intergranular  embrittlement  in  metals  is  usually 
caused  by  impurities  segregating  to  the  grain 
boundaries.  Impurities  present  in  bulk  concen¬ 
trations  of  only  10-100  parts  per  million  (ppm) 
can  result  in  a  dramatic  decrease  in  the  mech¬ 
anical  properties  (primarily  ductility,  fracture 
toughness  and  fracture  strength)  of  metallic 
materials,  and  thus  pose  significant  processing 
and  application  problems.  This  detrimental 
effect  of  mere  ppm  amounts  of  impurities  is 
consistent  with  the  fact  that  only  few  ppms  of 
impurities  are  sufficient  to  saturate  all  the  grain 
boundaries  in  a  polycrystal  of  the  typical  grain- 
size  (10-100  ^m)3.  The  sensitivity  of  DBTT  to 
the  grain  size  confirms  the  aforementioned 
effect  of  grain  boundary  impurities.  That  is,  the 
smaller  the  grain  size,  the  larger  the  amount  of 
impurities  required  to  saturate  the  grain  bound¬ 
aries  at  a  given  temperature  and  hence,  at  the 
same  bulk  concentration  of  impurities,  the 
higher  toughness  levels  and  the  lower  DB  — 
TT  are  expected  in  fine-grained  polycrystals. 

In  addition  to  the  smaller  grain  size,  increas¬ 
ing  the  purity  of  the  metallic  system,  in  general, 
increases  the  material’s  ductility  and  lowers  its 
DBTT.  For  example,  the  DBTT  in  high  purity 
W  single  crystals  processed  by  electron-beam 
zone  remelting,  involving  the  use  of  a  special 
impurity  gettering  procedure,  was  found  to  be 
as  low  as  —  196°C.  If  impurities  are  the  main 
cause  of  embrittlement  in  tungsten,  gettering 
the  impurities  is  the  obvious  way  of  ductilizing 
this  metal.  A  well-known,  though  extremely 
expensive,  purity  enhanced  procedure  is  the  so- 
called  ‘rhenium  effect’7  ”.  While  the  actual 
mechanism  of  the  rhenium  effect  is  not  yet 
completely  understood,  it  is  believed  to  be 
based  on  the  gettering  of  oxygen7,8.  A  more 


promising  way  of  removing  ‘harmful’  impurities, 
such  as  O,  N,  P,  etc.,  from  the  grain  boundary  is 
by  gettering  them  with  micro-alloying  additions 
of  selected  elements,  e.g.  boron,  which  leads  to 
the  formation  of  compounds  such  as  boron 
oxides,  boron  nitrides,  etc.g_u.  This  process, 
however,  requires  a  careful  control  of  the 
material  chemistry  since  ductility  improvements 
due  to  impurity  gettering  can  be  expected  only 
so  far  as  the  resulting  compound  precipitates 
are  fine.  The  coagulation  and  coarsening  of  the 
precipitates  generally  result  in  an  adverse 
embrittling  effect. 

using  both  semi-empirical  and  first-principles 
calculations,  Krasko3  carried  out  a  theoretical 
analysis  of  clean  tungsten  grain  boundaries  and 
grain  boundaries  containing  various  impurities. 
His  work  provided  some  important  results 
regarding  the  effect  of  various  grain  boundary 
impurities  and  micro-alloying  additions.  For 
instance,  impurities  such  as  H,  N,  O,  P,  S,  Si  are 
found  to  weaken  the  intergranular  cohesion  in 
tungsten.  On  the  contrary,  the  presence  of  B 
and  C  was  found  to  enhance  the  bonding  across 
the  grain  boundary,  thus  improving  the  inter¬ 
granular  cohesion.  Furthermore,  the  so-called 
site-competition  effect,  as  a  result  of  which  the 
species  with  a  lower  energy  at  the  grain  bound¬ 
ary  tend  to  replace  the  species  which  have  a 
higher  energy  at  the  grain  boundary,  was  found 
to  play  an  important  role  in  affecting  the 
impurity  distribution  at  the  tungsten  grain 
boundaries3.  Among  the  atomic  species  ana¬ 
lyzed,  B  was  found  to  have  the  lowest  energy  at 
the  grain  boundary,  and  thus  wound  tend  to 
displace  other  impurity  atoms  from  the  grain 
boundary,  while  at  the  same  time  enhancing  the 
intergranular  cohesion.  Based  on  these  findings, 
micro-alloying  of  technically  pure  tungsten  with 
10-50  ppm  of  B  was  recommended  as  a  way  of 
enhancing  the  ductility  of  this  metal3.  Experi¬ 
mental  investigation  of  tungsten  alloyed  with  B 
in  amounts  15  times  higher  than  the  one  recom¬ 
mended  by  Krasko3  was  found  to  result  in  a 
significant  (150°C)  drop  in  the  DBTT6-11. 
Microstructural  analysis  further  revealed  the 
presence  of  relatively  coarse  boron  oxide  par¬ 
ticles  which  are  the  result  of  the  excessive 
amount  of  B  in  the  tungsten.  This  suggests  that 
the  addition  of  10-50  ppm  of  B  which  is  suffi¬ 
cient  to  completely  displace  oxygen  from  the 
grain  boundaries  via  the  ‘site  competition’  effect 
without  an  excessive  formation  of  boron  oxide, 
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as  suggested  by  Krasko,  may  result  in  even 
higher  levels  of  ductility  and  fracture  toughness. 

In  the  present  work,  the  molecular  dynamics 
method  was  used  to  explore  the  effect  of  impur¬ 
ities  and  micro-alloying  additions  (B,  C,  N,  O, 
Al,  S,  Si,  P)  on  the  fracture  resistance  of  the 
S3  (111)  grain  boundary  in  tungsten.  The 
organization  of  the  paper  is  as  follows.  Modifi¬ 
cations  of  the  ‘environment  sensitive  embed¬ 
ding’  (ESE)  energy  functions  for  various  grain 
boundary  impurities  in  tungsten  are  discussed 
below.  Generation  of  the  computational  crystal 
containing  a  crack  along  the  S3  (111)  grain 
boundary  is  presented  in  the  section  ‘Computa¬ 
tional  crystal’.  Computations  of  the  ideal  work 
of  grain  boundary  separation  and  the  Fx  con¬ 
servation  integral  are  discussed  in  the  next  two 
sections,  respectively.  The  section  ‘Results  and 
Discussion’  contains  the  results  obtained  in  the 
present  study,  while  the  main  conclusions  are 
drawn  in  the  final  section. 

PROCEDURE 

Modification  of  the  Finnis-Sinclair  functions 
for  tungsten 

Within  the  EAM  formalism'2, 13  the  total  poten¬ 
tial  energy  of  the  system  is  given  as  the  sum  of 
two  energy  contributions,  the  interaction  of 
each  atom  with  the  local  electron  density  associ¬ 
ated  with  the  remaining  atoms  in  the  system, 
called  the  embedded  energy,  and  a  pair-wise 
interaction  term  reflecting  the  electrostatic 
interactions  between  the  atoms.  In  particular, 
the  total  potential  energy  is  written  as: 

£lot  =  LE,.(p,.)+^  E  <M*o)  (!) 

where,  F,  is  the  embedded  energy  of  atom  i,  p , 
is  the  electron  density  at  atom  i,  and  f/j,y(/?,y)  is 
the  pair-wise  interaction  between  atoms  i  and  j 
separated  by  the  distance  i?iy.  The  electron 
density  at  each  site  is  computed  from  the  super¬ 
position  of  spherically  averaged  atomic  electron 
densities  as: 

P,=  I  (2) 

pf(Ri})  in  eqn  (2)  represents  the  atomic  electron 
density  at  site  i  due  to  an  atom  at  site  j  at  a 
distance  Rir  The  superscript  a  in  p/(Riy)  is  used 
to  specify  the  species  of  the  atom  at  site  j.  The 


embedding  energies  and  pair-wise  interaction 
functions  in  eqn  (1)  are  generally  determined 
by  fitting  various  physical  quantities  of  the 
material  such  as  the  sublimation  energy,  the 
equilibrium  lattice  parameters,  the  second-order 
elastic  constants,  the  various  defect  formation 
energies  and  the  zero-temperature  equation  of 
state14. 

The  Finnis-Sinclair  method15  that  is  usually 
called  ‘the  N-body  potential  approach’  has  a 
slightly  different  foundation,  but  the  total 
potential  energy  can  be  still  expressed  in  the 
form  of  eqn  (1).  The  first  term,  however,  has  a 
different  meaning  since  it  originates  from  the 
tight-binding  theory  according  to  which: 

F  (pd  =  —  aJp,  (3) 

where  A  is  a  species-dependent  parameter 
whose  value  for  tungsten  (1-896373  eV)  was 
determined  by  Finnis  and  Sinclair15  and  p,  now 
represents  an  environment-dependent  effective 
atomic  coordination  number,  or  an  effective 
atomic  density  at  site  i.  Finnis  and  Sinclair15 
originally  proposed  the  following  functional 
form  for  p,: 

p'frR,,)  =  (RtJ - Fcu[)2,  for  R/j < Rcul 

p'frR.fr  =  0  for  R,,  >  /?cut  (4) 

where  the  superscript  h  stands  for  ‘host  atoms’. 
For  tungsten,  Finnis  and  Sinclair  set  the  cut-off 
distance  Rcut  to  a  value  4-400224  A.Finnis  and 
Sinclair15  further  proposed  the  following  form 
for  the  pair-wise  potential  function  </>,y: 

<l>u(Rn)  -  (Ru  -  cf(c0+clRu+c2Rfj),  for  R,,  <  c 
and 

PuiR.fr  =  0  for  Ru>c  (5) 

where,  for  tungsten,  the  cut-off  distance  c  is  set 
to  a  value  3-4  A,  and  the  quadratic  poly¬ 
nomial  coefficients  are  cG  =  47-13465,  cx  = 
—  33-7665655,  and  c2  =  6-254199.For  tungsten 
containing  interstitial  impurities,  Krasko3  intro¬ 
duced  a  modified  Finnis-Sinclair  formulation 
by  simply  adding  a  term  to  eqn  (1)  to  account 
for  the  energy  of  impurity  atoms  in  the  host 
environment.  The  modified  £tot  is  then  written 
as: 

Ftot  =  S  F(pD+^  X  <M*0-)  +  ?  ESEifrD  (6) 
'  2  k 
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where  first  two  terms  on  the  right-hand  side  of 
eqn  (6)  are  given  by  eqns  (3) -(5),  while  the 
third  term  represents  the  sum  of  the  ‘environ¬ 
ment-sensitive  embedding’  (ESE)  energies  of 
the  impurity  atoms,  and  is  equal  to  the  inter¬ 
action  energy  between  the  impurity  and  host 
atoms,  p*mp  is  the  electron  density  at  the  site  of 
impurity  atom  k  due  to  the  surrounding  host 
atoms,  and  is  given  as: 

prp=  i  pr=  e  pf(RkJ)  (7) 

j  j 

where  p,H(i?,y)  represents  the  electron  density 
defined  by  eqn  (4)  at  the  site  k  of  an  impurity 
atom  due  to  the  host  atom  at  site  j,  while  RkJ  is 
the  distances  between  the  impurity  atom  at  site 
k  and  a  host  atom  at  site  j. As  mentioned  above 
the  ESE  functions  introduced  by  Krasko3  were 
derived  in  the  spirit  of  EAM,  that  is  relative  to 
the  free  atom  electron  density  at  the  impurity 
site  due  to  the  host  atoms,  pj1,  while  the  host 
atom  embedded  function  are  expressed  relative 
to  the  effective  atomic  density,  p,h.  To  avoid 
possible  confusion  relative  to  the  meaning  of 
the  density  functions  used,  the  ESE  functions 
were  rederived  in  the  present  paper  relative  to 
the  effective  atomic  density  used  by  Finnis  and 
Sinclair,  p,h.  The  resulting  ESE(pimp)  expressed 
as  fourth-order  polynomials  are  given  as: 

ESE  (pimp)  =  p"""(C0+C1p,nip+C2(pimp)2+C3(pimp)3 

(8) 

The  values  of  the  coefficients  C0,  C,,  C2  and 
C3  for  various  impurities  and  micro- alloying 
additions  are  given  in  Table  1  and  the  corre¬ 
sponding  ESO  vs  pimp  functions  are  plotted  in 
Fig.  1. 


Table  1.  Environment-sensitive  embedded  energy,  ESE,  vs 
effective  atomic  density,  p,  function  coefficients  for 
various  impurities  in  tungsten 


ESE(eV) 

=  p(C  o+Crf+C^+C  3p3) 

Impurity  CQ 

c, 

c2 

c3 

B 

-0-8926 

0-0099 

-2-548x10  5 

7-833  xlO”9 

C 

-0-6177 

0-0125 

—  6-223  x  10  5 

1-124  xlO'7 

N 

-0-9423 

0-0157 

-6-222x10- 5 

7-834x10- 8 

O 

-0-6202 

0-0134 

-4-798x10- 5 

5-044x10- 8 

A1 

-0-6259 

0-0244 

-8-997  xlO5 

9-926x10  8 

Si 

-0-5201 

0-0233 

-8-199  xlO-5 

8-131  x  10  8 

S 

-0-1674 

0-0224 

-8-132  xlO-5 

8-522  xlO-8 

P 

-0-3655 

0-0299 

- 1-292  x  10  -4 

1-792x10- 7 

Fig.  1.  Environment  sensitive  embedding  (ESE)  energy  vs 
atomic  electron  density  functions  for  various  grain-bound- 
ary  impurities  and  micro-alloying  additions  in  tungsten. 


Computational  crystal 

Generation  of  the  grain  boundary 
simulations  were  carried  out  using  a  cylindrical 
crystal  whose  crystallographic  orientation  is 
shown  in  Fig.  2.  The  dimensions  of  the  crystal 
are  given  in  terms  of  the  number  of  (112), 
(111),  (110)  B.C.C.  interplanar  spacings  along 
the  three  principal  axes.  The  initial  diameter  of 
the  crystal  was  approximately  9  nm. 

To  generate  a  £3  (111)  grain  boundary,  the 
atoms  in  the  upper  part  of  the  computational 
crystal  (x2  >  0)  were  rearranged  to  produce  a 
configuration  which  is  a  mirror  image  of  the 


Fig.  2.  Schematic,  size  and  the  crystallographic  orientation 
of  the  computational  crystal  used  in  the  present  work. 
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lower  part  of  the  crystal  across  the  grain  bound¬ 
ary  plane  (x2  =  0).  The  resulting  structure  of  the 
computational  crystal  can  be  represented  by 
succession  of  the  (111)  B.C.C.  planes  as: 

CBACBACBABCABCABC  (8) 

where  the  grain  boundary  is  marked  by  A.  This 
grain  boundary  is  of  a  tilt  type  and  is  character¬ 
ized  by  an  inverse  coincidence  lattice  number 
Z  =  3.  Krasko3  showed  that,  due  to  the  associ¬ 
ated  lowest  energy  level,  an  impurity  atom  such 
as  B,  C,  N,  etc.,  is  most  likely  to  occupy  an 
interstitial  position  in  the  center  of  the  trigonal 
prism  formed  by  6W  atoms  in  the  Z3  (111) 
grain  boundary.  Consequently,  all  the  simula¬ 
tions  of  the  Z2  (111)  grain  boundary  containing 
impurities  were  done  under  the  condition  that 
the  impurity  atoms  occupy  the  interstitial  sites 
associated  with  the  trigonal  prisms  along  the 
grain  boundary  plane.  It  should  be  noted  that 
for  a  typical  grain  size  range  of  10-100  /im,  the 
introduction  of  interstitials  into  every  grain 
boundary  interstitial  site  corresponds  to  the 
impurity  concentration  range  of  80-20  ppm. 

Grain  boundary  crack 

To  generate  a  crack  along  the  S3  (111)  grain 
boundary,  all  the  atoms  in  the  computational 
crystal  were  displaced  from  their  initial  posi¬ 
tions  in  accordance  with  the  plane  strain  linear 
elastic  solution  developed  by  Sih  and  Liebo- 
witz16  for  the  crack  displacements  in  an  aniso¬ 
tropic  single  crystal.  While,  in  general,  the 
single-crystal  solution  of  Sih  and  Liebowitz16  is 
not  valid  for  interfacial  cracks,  the  fact  that  the 
two  crystals  in  Fig.  2  had  identical  correspond¬ 
ing  crystallographic  directions  in  the  three 
principal  directions,  justifies  the  use  of  this 
plane  procedure.  Under  a  pure  tensile  load 
applied  in  the  x2  direction  and  for  plane  strain 
condition  along  the  x3  direction  the  components 
of  the  displacement  along  the  three  principal 
axes  are  given  as  following: 


u,=Kt 


2  r 

—  Re 

7t 


1 


S,  S2 


x  [j,p2(cos  9+s2  sin  9)1' 2 


—  s2jp,(cos  9+s,  sin  dy 


w2  —  K  j 


2  r 

—  Re 
n 


1 


S i  S2 


x  [s,42(cos  9+s2  sin  0),/2 


—  s2qt( cos  9+s{  sin  0)'/2] 


«3  =  0  (9) 

where  r  and  9  are  the  polar  coordinates,  i.e.  the 
radial  distance  and  the  angle  between  the  radial 
line  and  the  Xj-axis  in  Fig.  2.  Functions  p,  and 
q}{j  =  1,2)  are  defined  as: 

Pi  =  tti  i s j ~t~a 1 2  Oi&S] 

p2  —  a !  iS2+a]2  al6s2 

at2Si  +a22  a26s  t 
<h= - 

Si 

ai2s2+a22  a2es2 

42=  -  (10) 

s2 


where  s,  =  +//?,  and  s2  =  fi2  =  ot2+i(}2  are 

two  non-complex  conjugate  roots  of  the  follow¬ 
ing  characteristic  equation: 

a , ,  p*  -  2a ,  6ju3+(2a ,  2+«66)^2  -  2 a26pj+a22  =  0  (11) 


The  remaining  two  solutions  of  eqn  (11)  are  the 
corresponding  complex  conjugates  of  the  first 
two,  ie  p3  =  n*1  and  f4  =  f*2.  coefficients 
appearing  in  eqns  (10)  and  (11)  are  related  to 
the  elastic  compliance  constants  of  the  material, 
Sij,  and  for  plane-strain  are  given  as: 


a22  — 


‘^22^33  *^23 


^11^33  *^13  ^12^33  *^13^3 

O’  11=  “  >  U  1 2  =  fl2 1  =  ~ 

O33  J33 

Sl6$33  ~Sl3$36 

S33 

$26$  33  $  1 3^36 

S33 


"  *  Cl  16  —  Clf,  1  — 


U(s<s  — 


$66$  33  $36 


"  1  a26  —  ae2  : 


(12) 

For  a  given  material,  the  elastic  compliance 
constants  used  in  eqn  (12)  depend  on  the 
crystallographic  orientation  of  the  computa¬ 
tional  crystal  being  examined.  Parameter  K, 
appearing  in  eqn  (9)  is  the  mode  I  stress  inten¬ 
sity  factor  whose  critical  value  associated  with 
the  reversible  crack  extension,  the  Griffith  stress 
intensity  factor,  KGi,  is  given  as16: 
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Kar  =  j2ysep/A  (13) 

where  A  is  defined  as: 


A  = 


&\  1^22  \ 

/  a22 

\  ,  /2  2fl  i  2+^66 

2  J 

\  a, , 

/  2a, ,  j 

(14) 


and  this  parameter  also  in  the  relationship 
between  the  /,-integral  and  the  stress  intensity 
factor  Kj,  Jl=AK$  [eg  16].The  2ysep  term 
appearing  in  eqn  (13)  represents  the  ideal  work 
of  grain  boundary  separation  and  is  equal  to  the 
amount  by  which  the  work  per  unit  area  of  the 
crack  surface  done  by  the  external  loads 
exceeds  the  change  in  the  elastic  strain  energy. 
At  the  onset  of  crack  growth,  2ysep  is  given  as: 

2y«P  =  tf+tf-yZ,B  (15) 

where  superscripts  A  and  B  denote  the  two 
crystals  being  jointed  along  the  grain  boundary, 
ys  is  the  surface  energy  and  yint  the  grain  bound¬ 
ary  (interface)  energy. 


Computation  of  the  ideal  work  of  grain 
boundary  separation,  2yint 


To  compute  the  surface  energy  associated  with 
the  (lll)bcc  crystallographic  plane,  the  crack 
plane,  a  rectangular  slab  with  the  same  orienta¬ 
tion  as  the  computational  crystal  shown  in  Fig. 
2  and  with  its  top  and  bottom  surfaces  parallel 
to  this  plane  is  created.  The  height  of  the 
crystal  is  chosen  to  be  large  enough  to  avoid  the 
interaction  between  the  top  and  the  bottom  sur¬ 
faces.  With  the  periodic  boundary  conditions 
applied  in  the  other  two  orthogonal  directions, 
the  energy  of  the  slab  is  minimized.  The  surface 
energy  is  then  defined  as  the  excess  energy  per 
unit  surface  area  of  this  crystal,  relative  to  the 
one  of  the  same  size  in  which  the  periodic 
boundary  conditions  are  applied  in  all  three 
directions. 

The  above  procedure  is  repeated  using  a  slab 
containing  the  S3  (111)  grain  boundary  and  the 
corresponding  grain  boundary  energy  calculated 
as  the  excess  energy  of  the  slab  per  unit  grain 
boundary  area  relative  to  the  perfect  B.C.C. 
without  a  grain  boundary  subject  to  the  periodic 
boundary  condition  in  all  the  direction. 

The  same  procedure  for  calculating  the  sur¬ 
face  and  the  grain  boundary  energy  was  utilized 
for  the  cases  where  the  initial  grain  boundary 


Table  2.  (lll)bcc  surface  energy,  ysur(,  £3  (lll)bcc  grain 
boundary  energy,  }>GB,  and  the  corresponding  ideal  work 
of  grain  boundary  decohesion,  2yint,  for  tungsten  contain¬ 
ing  various  impurities 


Impurity 

7surg,  (eV/A2) 

7ob,  (eV/A2) 

2yint,  (eV/A2) 

pureW 

0-4285 

0-5823 

0-2747 

B 

0-3006 

0-0909 

0-6382 

C 

0-3007 

0-1216 

0-6076 

N 

0-2583 

0-0943 

0-5925 

O 

0-2135 

0-3290 

0-3130 

A1 

0-1746 

0-3263 

0-2768 

Si 

0-1925 

0-3074 

0-3136 

S 

0-0735 

0-3515 

0-1505 

P 

0-1100 

0-3275 

0-2110 

contains  one  of  the  impurities  or  micro-alloying 
additions.  The  results  of  the  above  procedure 
for  pure  tungsten  and  tungsten  containing 
impurities  are  summarized  in  Table  2. 

Molecular  dynamics  method 

The  evolution  of  the  material  surrounding  the 
crack  tip  is  studied  using  the  standard  molecu¬ 
lar  dynamics  procedure  at  100  K.  The  molecular 
dynamics  calculations  involve  the  solution  of  the 
Newton  equations  of  motion  for  a  system  of 
interacting  particles  (atoms)  in  order  to  deter¬ 
mine  the  classical  particles  trajectories  and  velo¬ 
cities. 

At  the  beginning  of  each  simulation  run,  a 
crack  was  generated  by  displacing  all  atoms 
according  to  eqn  (9)  from  their  initial  positions 
corresponding  to  the  perfect  unrelaxed  B.C.C. 
bicrystal  structure  shown  in  Fig.  2.  All  the 
molecular  dynamics  simulation  runs  were  done 
under  the  fixed  periodic  boundary  conditions  in 
the  x3  =  [110]bcc  direction  to  ensure  the  plane 
strain  condition  in  this  direction.  Also,  at  each 
level  of  Ku  fixed  displacement  boundary  condi¬ 
tions  were  applied  onto  the  outer  surface  of  the 
computational  crystal. 

The  temperature  at  which  all  the  simulation 
runs  were  carried  out,  100  K,  was  set  by  initially 
assigning  to  each  atom  a  random  velocity 
according  to  the  Boltzmann  distribution.  During 
the  simulatino  process,  the  temperature  was 
maintained  constant  by  exponentially  relaxing, 
at  each  time  step  (2  fs),  the  average  squared 
atomic  velocity  using  a  time  constant  of  0T  ps. 
Using  this  procedure,  the  temperature  could  be 
maintained  within  3%  of  its  target  value. 
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Computation  of  F,  integral 

To  quantify  the  effect  of  the  crack-tip  stress 
relaxation  processes  on  the  resistance  of  a  S3 
(111)  grain  boundary  to  separation,  the  Eshel- 
by’s  conservation  F  integral17  was  calculated  for 
both  the  pure  tungsten  case  and  the  cases  of 
tungsten  containing  one  of  the  impurities  or 
micro-alloying  additions.  The  F  integral  pro¬ 
vides  a  means  for  determination  of  the  energy 
release  rate  accompanying  the  crack  extension 
in  cases  where  material  nonlinearity  effects  can¬ 
not  be  neglected.  F  is  a  vector  and  its  compo¬ 
nents  along  the  three  principal  axes  (the 
cracking  direction  xu  the  crack  plane  normal 
direction  x2  and  the  crack  front  direction  x3 ) 
represent  the  three  force  components  acting  on 
the  crack  tip.  For  the  computational  crystal  and 
the  crack  orientation  used  in  the  present  work 
(Fig.  2),  F,  acts  to  propagate  the  crack  tip, 
while  F2  and  F3  are  not  physical,  the  compo¬ 
nents  of  F  can  be  defined  using  a  closed  con¬ 
tour  T  which  surrounds  the  crack  tip  as 
following:17 


F,= 


kj 


Quk 

0  Xi 


d  Sj 


(16) 


as  being  quadratic  in  incremental  strains,  the 
stresses  are  calculated  using  the  following  cen¬ 
tral  difference  formula: 


dW  WiXf+dd-Wfa-S,) 
3e(/  2  A  u,j 


(18) 


where  etj  are  the  components  of  the  strain  ten¬ 
sor,  <5,  =  u^j  is  a  small  perturbation  of  the 
atomic  position  which  is  the  result  of  the  appli¬ 
cation  of  a  small  uniform  displacement  gradient 
Ally  to  the  entire  crystal.  The  value  of 
A Uij  =  0-0005  was  found  to  be  optimal  in  the 
present  work.  To  compute  the  components  of 
the  stress,  one  of  the  A utj  is  set  to  0-0005  at  a 
time  (the  other  Su are  kept  at  zero),  and  corre¬ 
sponding  nonzero  <5,  computed  and  used  in  eqn 
(18). 

From  the  theory  of  elasticity  it  is  known  that 
the  distortions  w,y  =  dujdxj  are  composed  of 
strain  (e,y)  and  rotational  (w/;)  components,  i.e. 

=  Ey+Wy.  A  strain  ellipsoid  that  relates  the 
extensional  strain  between  an  atom  and  its  a-th 
neighbor:  ea  =  (r“)  —  1,  where  rJ  and  R  are 
respectively  the  atomic  distances  after  and 
before  straining,  can  be  used  to  obtain  the  local 
strain  components  e(/  as  following: 


:  =  /°7*e.. 

li  lj*zij 


(19) 


where  W  is  the  strain  energy  density,  uk(xf)  is 
the  k-component  ( k  =  1,2,3)  of  the  displacement 
at  a  point  represented  by  the  coordinates 
Xi(i  =  1,2,3),  <r,y  are  the  stress  components  and 
dSj  =  dS-n,  where  is  the  j  component 
(/  =  1,2,3)  of  the  unit  outward  normal  vector  to 
the  contour  segment  of  length  dS.To  calculate 
the  stresses,  the  procedure  proposed  by  Hoag- 
land  et  al.18  was  used.  According  to  this  pro¬ 
cedure  the  strain  energy  density  associated  with 
an  atom  at  the  position  x,  is  defined  as: 


H/(x,)  = 


E(Xi)-E0 


Q 


O 


(17) 


where  E(xf)  is  the  potential  energy  of  the  atom 
with  the  coordinates  x,(«  =  1,2,3),  E0  is  the  equi¬ 
librium  energy  which  the  atom  would  have  in  a 
perfect,  stress-free  bulk  crystal  and  Q0  is  the 
equilibrium  atomic  volume.  The  energies  E(xf) 
and  E0  are  calculated  using  the  Finnis-Sinclair 
interatomic  potentials  and  an  energy  expression 
for  the  potential  per  atom  analogous  to  eqn 
(1). Assuming  that  for  small  strain  increments 
the  strain  energy  density  can  be  approximated 


where  If  are  the  directional  cosines  (in  the 
unstrained  lattice)  to  a  atom,  eqn  (19)  simply 
represents  a  transformation  of  the  strain  com¬ 
ponents  from  the  coordinate  system  given  in 
Fig.  2,  to  the  one  whose  xx  axis  is  along  the  line 
connecting  the  two  atoms  in  question,  i.e.  the 
atom  at  whose  location  the  strain  components 
are  being  evaluated  and  its  a-neighbor.  To  find 
the  strain  components  e,,  at  the  position  of  an 
atom  at  least-squares  procedure  was  used  to 
minimize  the  following  sum  over  the  n  neigh¬ 
bors  of  each  atom: 

Le=  £  (e“-/“/*e,v)2  (20) 

a  —  1 

The  minimization  procedure  dLJdSij  =  0,  yields 
six  linear  algebraic  (normal)  equations  which 
were  readily  solved  using  the  gauss  elimination 
method  (e.g.  19).  The  basis  for  determination  of 
rigid  body  rotations,  wjy,  at  a  particular  site  is 
obtained  by  proceeding  in  a  manner  similar  to 
the  one  employed  for  the  strains.  The  relation¬ 
ship  between  displacements  of  the  nearest 
neighbors  from  the  given  atom,  uf,  and  the  anti- 
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symmetric  matrix  of  the  local  rotations  suggests 
that  the  following  sum  should  be  minimized: 

LR  =  z  [(«T+w ,  2*2  -  w ,  3xlf+(ul  -  w ,  2x*+w23XiY 

a 

+  («3+w ,  3x*  —  W23X2)2]  (21) 

The  minimization  procedure  dL^dw^  =  0  yields 
three  normal  equations  which  were  also  readily 
solved  using  the  gauss  elimination  method.The 
evaluation  of  the  F,  integral  given  by  eqn  (16) 
was  done  numerically  in  the  present  work  using 
the  trapezoidal  rule19  along  a  circular  contour 
around  the  crack.  Owing  to  the  discrete  nature 
of  the  crystal,  all  the  atoms  within  a  distance  of 
0-2  nm  from  the  contour  were  assumed  to  be 
associated  with  the  contour.  The  magnitude  of 
the  contour  segment  dSa  corresponding  to  an 
atom  a  is  calculated  as  dS“  =  7?(0a+1  —  0a+1)/2 
where  R  is  the  contour  radius,  and  6a+1  and 
0,  _  1  are,  respectively,  the  polar  angles  for  the 
atoms  in  the  contour  which  immediately  follow 
and  proceed  atom  a  in  the  counterclockwise 
direction.  The  unit  vector  normal  n01  to  the  con¬ 
tour  segment  d 5a  is  defined  as  =  cos  6J- 
-l-sin  6J,  where  i  and  j  are  the  unit  vectors  in  x1 
and  x2  directions,  respectively.20 

It  should  be  noted  that  the  strain  energy 
density  W,  and  hence  F,  each  have  two  con¬ 
tributions:  one  associated  with  the  strain  energy 
density  and  the  other  arising  from  the  work  of 
grain  boundary  separation.  The  former  con¬ 
tribution  to  F1;  and  /,  integral,  acts  to  extend 
while  the  latter,  2ysep,  acts  to  close  the  crack. 
Hence,  the  following  relationship  holds: 

Fi  =  7 ,  —  2ysep  (22) 

RESULTS  AND  DISCUSSION 

Ideal  work  of  grain  boundary  decohesion 

The  results  of  our  calculation  of  the  (111)  sur¬ 
face  energy,  E3  (111)  grain  boundary  energy 
and  of  the  corresponding  work  of  grain  bound¬ 
ary  decohesion  are  listed  in  Table  2.  The  results 
pertain  to  the  case  when,  after  decohesion,  all 
the  impurity  atoms  remain  on  one  of  the  crack 
faces  whereas  the  other  crack  face  is  clean  of 
impurities.  While  the  magnitudes  of  the  work  of 
separation  takes  on  different  values  when  the 
impurities  are  assumed  to  divide  equally 
between  the  two  crack  surfaces  (data  not 
shown),  the  general  findings  regarding  the 


effect  of  various  grain  boundary  impurities  on 
the  grain  boundary  decohesion  resistance 
remain.  The  results  shown  in  Table  2  indicate 
that,  on  the  one  hand,  elements  such  as  B,  C 
and  N  increase  the  work  of  E3  (111)  grain 
boundary  decohesion.  On  the  other  hand,  ele¬ 
ments  such  as  O,  A1  and  Si  have  a  relatively 
small  effect  on  this  work.  Lastly,  elements  such 
as  P  and  S  act  as  strong  grain  boundary  embrit- 
tlers  and  reduce  this  work.  These  findings  are 
further  correlated  with  the  grain  boundary 
energy  data  shown  in  Table  2.  The  grain  bound¬ 
ary  energy  data  shown  that  B,  C  and  N  reduce 
the  grain  boundary  energy  the  most,  while  P 
and  S  reduce  this  energy  the  least.  This  suggests 
that  B,  C  and  N  would,  via  the  site-competition 
effect,  tend  to  displace  P  and  S  from  the  grain 
boundary,  and  in  turn,  enhance  the  grain 
boundary  cohesion.  In  other  words,  the  deleteri¬ 
ous  embrittling  effect  of  impurities  such  as  P 
and  S  in  tungsten  can  be,  as  previously  sug¬ 
gested  by  Krasko,3  reduced  by  using  micro¬ 
alloying  additions  of  micro-alloying  elements 
such  as  B,  C  and  N. 

Initial  atomic  configurations 

The  initial  atomic  configuration  of  the  Z3  (111) 
grain  boundary  in  pure  tungsten  (‘clean’  grain 
boundary)  containing  a  crack  at  the  applied 
stress  intensity  factor  KHpp  =  2-0  KGt  =  2-7694 
MPa  m1/2  and  the  corresponding  configurations 
after  relaxation  (energy  minimization)  are 
shown  in  Fig.  3(a)  and  (b),  respectively.  The 
energy  minimization  was  carried  out  using  the 
conjugate  gradient  method.21  The  structure  of 
the  grain  boundary  ahead  of  the  crack  after 
relaxation  can  be  characterized  by  analyzing  the 
spacing  of  the  adjacent  (111)  planes  as  a  func¬ 
tion  of  the  distance  from  the  grain  boundary 
and  comparing  it  with  the  bulk  value  in  pure 
tungsten  d(111)  =  0-9137  A.  As  shown  in  Table  3, 
the  (111)  interplanar  spacing  oscillates  with  the 
distance  from  the  grain  boundary  with  the 
amplitude  of  oscillations  gradually  decreasing. 
Consequently  by  the  10- 12th  plane  away  from 
the  grain  boundary  (not  shown  in  Table  3),  the 
oscillations  are  practically  damped  out.  The 
results  shown  in  Table  3  correspond  to  the  aver¬ 
age  (111)  interplanar  spacings  at  a  distance 
between  16-9  and  19-6  A  from  the  crack  tip  in 
the  positive  x1  direction.  It  should  be  noted  that 
the  distance  between  the  second  and  the  third 
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(a)  (b) 

Fig.  3.  Initial  atomic  configuration  of  the  ‘clean’  13  (111)  grain  boundary  in  tungsten  containing  a  crack:  (a)  before 
relaxation,  and  (b)  after  relaxation.  The  applied  level  of  the  stress  intensity  factor  Kapp  =  2-0  Kar  =  2-7694  MPa  m1/2. 


Table  3.  The  spacing,  in  A,  between  two  adjacent  (111) 
planes  as  a  function  of  the  distance  from  the  grain  bound¬ 
ary  in  pure  tungsten  and  in  tungsten  containing  B  and  P 
grain  boundary  impurities 


Adjacent  (111) 
planes 

Type  of  grain  boundary 

Clean 

With  B 

With  P 

Above  the  grain  boundary  plane 
G.B.  Plane-Plane  1  1-292 

1-334 

1-555 

Plane  1 -Plane  2 

0-609 

0-801 

0-906 

Plane  2-Plane  3 

1-009 

0-828 

0-734 

Plane  3-Plane  4 

0-984 

0-987 

0-977 

Plane  4-Plane  5 

0-824 

0-865 

0-906 

Below  the  grain  boundary  plane 

G.B.  Plane-Plane  1  1-252 

1-338 

1-570 

Plane  1-Plane  2 

0-566 

0-808 

0-908 

Plane  2-Plane  3 

1-051 

0-830 

0-734 

Plane  3-Plane  4 

0-973 

0-987 

0-968 

Plane  4-Plane  5 

0-821 

0-863 

0-917 

plane  has  decreased  significantly  relative  to  that 
in  the  bulk  B.C.C.  tungsten.  This  behavior  is  the 
manifestation  of  the  crystal’s  tendency  to 
reduce  the  free-volume  of  the  grain  boundary. 
It  should  also  be  noted  that  because  the  plane 
of  the  crack  lies  between  the  grain  boundary 
plane  and  the  first  (111)  plane  below  the  grain 
boundary  plane,  the  variation  of  the  (111)  inter- 
planar  spacings  with  the  distance  from  the  grain 
boundary  is  different  for  the  (111)  planes  below 
and  those  above  the  grain  boundary  plane. 

The  initial  atomic  configuration  of  the  S3 
(111)  grain  boundary  with  a  crack  in  tungsten 
containing  boron  impurity  atoms  located  in  trig¬ 
onal  interstitial  sites  and  the  corresponding  con¬ 
figuration  after  energy  minimization  are  shown 
in  Fig.  4(a)  and  (b),  respectively.  The  relaxed 
configurations  for  the  same  type  of  grain 


(a)  (b) 

Fig.  4.  Initial  atomic  configuration  of  the  S3  (111)  grian  boundary  in  tungsten  containing  a  crack  with  B  atoms  (smaller 
balls)  located  in  trigonal  interstitial  grain  boundary  sites:  (a)  before  relaxation,  and  (b)  after  relaxation.  The  applied  level 

of  stress  intensity  factor  Kapp  =  1-4  KCi,  =  2-9548  MPa  m1/z. 
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boundary  containing  other  grain  boundary 
cohesion  enhancers  such  as  C  and  N  are  quite 
similar  to  the  one  shown  in  Fig.  4(b)  and  are 
not  shown  here  for  brevity.  The  variation  of  the 
(111)  interplanar  spacings  with  distance  from 
the  grain  boundary  for  the  relaxed  configuration 
depicted  in  Fig.  4(b)  is  listed  in  Table  3.  In 
comparison  with  the  ‘clean’  grain  boundary 
beyond  the  two  planes  adjacent  to  the  grain 
boundary  plane,  the  magnitudes  of  the  oscilla¬ 
tion  amplitude  of  the  (111)  interplanar  spacings 
is  significantly  lower.  In  addition,  the  distance 
of  the  two  nearest  (111)  and  the  two  next-near¬ 
est  (111)  planes  to  the  grain  boundary  increase 
somewhat  due  to  the  presence  of  interstitial 
impurities  of  B,  C  or  N  relative  to  those  in  the 
case  of  clean  grain  boundary  in  tungsten,  Table 
3. 

The  initial  atomic  configuration  of  the  S3 
(111)  grain  boundary  with  a  crack  in  tungsten 
containing  phosphorus  impurity  atoms  located 
in  trigonal  interstitial  sites  and  the  correspond¬ 
ing  configuration  after  energy  minimization  are 
shown  in  Fig.  5(a)  and  (b),  respectively.  The 
relaxed  configuration  for  the  same  type  of  grain 
boundary  containing  other  grain  boundary 
embrittlers  such  as  S  is  quite  similar  to  the  one 
shown  in  Fig.  5(b)  and  is  not  shown  here  for 
brevity.  The  variation  of  the  (111)  interplanar 
spacings  with  distance  from  the  grain  boundary 
for  the  relaxed  configuration  depicted  in  Fig. 
5(b)  is  listed  in  Table  3.  In  comparison  with  the 
‘clean’  grain  boundary,  the  magnitudes  of  the 
oscillation  amplitude  is  significantly  larger.  Fur¬ 
thermore,  the  distance  of  the  two  nearest  and 


the  two  next-nearest  (111)  planes  from  the 
grain  boundary  is  significantly  increased  relative 
to  those  observed  in  the  cases  of  the  clean  grain 
boundary  and  the  grain  boundary  containing  B, 
C  or  N  impurities. 

Material  evolution  at  the  crack  tip 

The  progress  of  material  evolution  at  the  crack 
tip  in  the  case  of  a  clean  S3  (111)  grain  bound¬ 
ary  at  the  applied  stress  intensity  factor 
^app  =  2-0  Kar  =  2-7694  MPa  m1/2  is  shown  in 
Fig.  6(a)  and  (b).  This  level  of  the  applied  K,  is 
the  minimal  level  at  which  a  clear  evidence  of 
crack  propagation  by  05  ps  of  the  molecular 
dynamics  simulation  time  can  be  obtained.  A 
comparison  of  the  atomic  configuration  after 
05  ps  simulation  time,  Fig.  6(a),  with  the  origi¬ 
nal  configuration,  Fig.  3(a),  shows  a  significant 
rearrangement  of  the  atoms  in  the  region  sur¬ 
rounding  the  crack  tip.  This  rearrangement 
leads  to  a  significant  crack  blunting.  No 
evidence  of  dislocation  emission  from  the  crack 
tip  can  be  found.  Atomic  rearrangement  with  a 
small  additional  crack  blunting  continues  to 
take  place  with  the  simulation  time  which  can 
be  readily  established  by  analyzing  the  atomic 
configuration  after  5  ps,  Fig.  6(b).  It  is  interest¬ 
ing  to  note  that  as  a  result  of  the  aforemen¬ 
tioned  material  evolution,  a  one-atom  thick, 
one-lattice  parameter-wide  strip  of  the  (110)bcc 
plane  has  formed  at  the  very  crack  tip,  Fig. 
6(b).  This  plane  is  normal  to  the  original  crack 
plane,  (111),  and  a  [100]  direction  in  this  plane 
is  aligned  with  the  initial  crack  front-direction 


(a)  (b) 


Fig.  5.  Initial  atomic  configuration  of  the  S3  (111)  grain  boundary  in  tungsten  containing  a  crack  with  P  atoms  (smaller 
balls)  located  in  trigonal  interstitial  grain  boundary  sites:  (a)  before  relaxation,  and  (b)  after  relaxation.  The  applied  level 

of  stress  intensity  factor  K„vp  =  1-2  Kar  =  T4563  MPa  m1/2. 
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[110] .  The  formation  of  this  strip  not  only 
causes  additional  crack  blunting,  but  also 
changes  the  local  ideal  work  of  grain  boundary 
decohesion,  which  in  this  case  can  be  approxi¬ 
mated  as  2y' '1<,l!,),  where  y^'1,1,0'  is  the  surface 
energy  of  the  (100)  plane.  Using  the  aforemen¬ 
tioned  procedure  for  the  calculation  of  the  sur¬ 
face  energy,  the  ideal  work  of  grain  boundary 
decohesiono  has  been  determined  as  2y^cr,f,>  = 
0-8128  eV/A2.  Since  this  value  is  higher  than  the 
one  0-2747  eV/A2  reported  previously  in  Table 
2,  this  may  be  the  reason  why  crack  propagation 
stops  once  the  (110)  strip  forms  at  the  grain 
boundary. 

The  material  evolution  in  the  region  sur¬ 
rounding  the  crack  tip  in  the  cases  of  the  S3 

(111)  grain  boundary  containing  interstitially 


dissolved  O,  A1  or  S  atoms  is  quite  similar  to 
that  of  the  clean  grain  boundary  and  will  not  be 
discussed  any  further. 

The  progress  of  material  evolution  at  the 
crack  tip  in  the  case  of  S3  (111)  grain  boundary 
containing  interstitially  dissolved  B  atoms  at 
the  applied  stress  intensity  factor  Kapp  = 
1-6  Kar  =  3-3770  MPa  m 1/2  is  shown  in  Fig.  7(a) 
and  (b).  K,  =  0-8  KGt  =  1-6885  MPa  m1/2  is  the 
minimal  Kt  value  at  which  a  noticeable  crack 
advancement  could  be  found  by  0-5  ps  of  the 
molecular  dynamics  simulation  time.  A  com¬ 
parison  of  the  atomic  configuration  after  0-5  ps 
simulation  time,  Fig.  7(a),  with  the  correspond¬ 
ing  original  configuration,  Fig.  4(a)  shows  that  a 
significant  material  evolution  takes  place  in  the 
region  surrounding  the  crack  tip.  This  evolution 
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gives  rise  to  a  major  crack  tip  blunting,  but  no 
evidence  of  dislocation  emission  can  be  found. 
The  aforementioned  material  evolution  and  the 
accompanying  crack-tip  blunting  continues  with 
the  simulation  time,  Fig.  7(b).  It  is  also  import¬ 
ant  to  note  that  after  5  ps,  Fig.  7(b),  B  atoms 
form  a  cluster  at  the  crack  tip.  This  cluster  is 
expected  to  enhance  the  intergrain  bonding 
and,  in  turn,  the  ideal  work  of  work  of  grain 
boundary  decohesion.  This  may  be  the  reason 
why  once  the  cluster  of  B  atoms  is  formed,  the 
crack  propagation  ceases. 

The  progress  of  material  evolution  at  the 
crack  tip  in  the  case  of  the  £3  (111)  grain 
boundary  containing  interstitially  dissolved  P 
atoms  at  the  applied  stress  intensity  factor 
Kapp  =  l'2Kar  =  1-4563  MPa  m1/2  is  shown  in 
Fig.  8(a)  and  (b).  K,  =  0-4  Kai  =  0-4854  M- 
Pa  m1/2  is  the  minimal  value  of  the  stress  inten¬ 
sity  factor  at  which  a  noticeable  crack  extension 
could  be  observed  by  0-5  ps  of  the  simulation 
time.  A  comparison  of  the  atomic  configuration 
after  0-5  ps  of  the  simulation  time,  Fig.  8(a), 
with  the  corresponding  initial  configuration  Fig. 
5(a),  shows  that  a  significant  rearrangement  of 
the  atoms  takes  place  in  the  region  of  the  grain 
boundary  ahead  of  the  crack.  This  evolution 
causes  a  significant  increase  in  the  interatomic 
distances  in  the  grain  boundary  region.  This, 
however,  does  not  cause  any  crack  tip  blunting, 
and  the  crack  tip  remains  as  sharp  as  in  the 
original  configuration.  The  aforementioned 
material  evolution  continues  with  the  simulation 
time  and  after  5  ps,  Fig.  8(b),  the  crack  tip  has 


advanced  by  approximately  17  A  relative  to  its 
original  position.  If  the  simulation  time  is 
increased  beyond  5  ps,  the  crack  continues  to 
advance  until  it  reaches  fixed  atoms  at  the  outer 
rim  of  the  computational  crystal. 

Mechanics  of  grain  boundary  fracture 

To  qualify  the  effect  of  the  various  material 
evolution  processes  described  in  the  previous 
section  on  the  intrinsic  resistance  of  the  £3 
(111)  grain  boundary  in  pure  tungsten  and  in 
tungsten  containing  different  impurities  micro¬ 
alloying  additions,  the  Eshelby’s  F,  conservation 
integral  was  calculated.  As  discussed  above,  F, 
represents  the  force  acting  on  the  crack  tip  in 
the  x1  direction.  The  time  evolution  of  this  force 
whose  initial  value  corresponds  to  the  minimal 
level  of  the  applied  stress  intensity  factor  at 
which  a  noticeable  crack  advance  could  be  seen 
by  0-5  ps  simulation  time  is  given  in  Table  4  and 
plotted  in  Fig.  9. 

In  each  of  the  cases  analyzed,  a  positive  value 
of  Fj  is  required  to  initiate  crack  propagation.  It 
appears  that  the  initial  values  of  F\  fall  into 
three  groups:  o  the  highest  values 

(F1«l-4±0-l  eV/A2)  for  the  £3  (111)  grain 
boundary  containing  interstitially  dissolved  B,  C 
or  N  atoms;  the  intermediate  value 
(F,  « 1  -0  +  0- 1  e V/A2)  for  the  clean  £3  (111) 
grain  boundary  and  the  one  containing  dis¬ 
solved  O,  A1  or  Si  atoms  and  the  lowest  value 
(Fj^O-7  +  OT  eV/A2)  for  the  grain  boundary 
containing  P  or  S  impurities.  These  findings  are 


(a)  (b) 

Fig.  8.  Atomic  configuration  of  the  E3  (111)  grain  boundary  in  tungsten  containing  phosphorus  atoms  (smaller  balls)  after 
the  molecular  dynamics  simulation  times  of:  (a)  0-5  ps,  and  (b)  5  ps.  The  applied  stress  intensity  factor 

Kapp  =  1-2  Kar  =  1-4563  MPa  mI/2. 
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Table  4.  The  stress  intensity  factor  for  a  reversible  extension  of  the  grain  boundary  crack,  Kar  [eqn  ()12],  the  critical  stress 
intensity  factor  for  crack  extension  obtained  during  molecular  dynamics  simulations,  Kc„t,  and  the  corresponding  value  of 
the  Eshelby’s  F,  conservation  integral 


Impurity 

KGn  ( MPa  m"2) 

K„it,  (MPa  m"2) 

0  ps 

F„ 

0-5  ps 

(eV/A2) 

1-0  ps 

pure  W 

1-3847 

1-6616 

0-8397 

B 

2-1106 

1-6885 

1-4431 

0-6313 

C 

1-8528 

0-6127 

N 

1-9320 

1-3836 

O 

1-4781 

1-4781 

0-4424 

A1 

0-4137 

Si 

1-4795 

0-5918 

0-3857 

S 

0-4100 

0-6973 

0-4944 

P 

1-2136 

0-4854 

0-8257 

0-6612 

0-5485 

fully  consistent  with  the  results  for  the  ideal 
work  of  grain  boundary  decohesion  listed  in 
Table  2. 

Fig.  9  further  shows  that  as  the  simulation 
time  increases  the  aforementioned  material 
evoat  the  crack  tip  continues  to  take  place  and, 
as  a  result,  the  magnitude  of  the  crack  tip 
decreases.  This  decrease  between  0  and  1  ps  is 
significant  in  all  the  cases  except  for  the  ones 
involving  interstitially  dissolved  P  and  S  atoms 
at  the  grain  boundary  where  the  decrease  in  Fx 
with  the  simulation  time  is  quite  small.  Past  1  ps 
no  major  change  in  the  magnitude  of  F,  is 
found.  In  other  words,  the  material  evolution  at 
the  crack  tip  causes  F,  to  decrease  but  only  to  a 


lower  positive  value.  In  the  cases  of  the  clean 
S3  (111)  grain  boundary  and  such  a  boundary 
containing  B,  C,  N,  O,  A1  or  Si,  this  lower  level 
of  Fj  is  not  sufficient  to  give  rise  to  any  addi¬ 
tional  crack  extension  and  consequently,  the 
crack  growth  ceases.  In  sharp  contrast,  when 
the  grain  boundary  contains  P  and  S  impurities 
the  crack  growth  continues  while  the  F,  remains 
fairly  constant.  The  continued  crack  growth  in 
the  case  of  the  grain  boundary  containing  P  and 
S  atoms  is  promoted  by  a  continued  rearrange¬ 
ment  of  the  grain  boundary  atoms  which  causes 
an  increase  in  the  interatomic  distances.  In 
other  words,  the  atomic  density  in  the  grain 
boundary  region  continues  to  decrease,  and  in 
turn,  lowers  the  cohesion  strength  of  the  grain 
boundary. 

In  summary,  the  atomic  simulation  results 
obtained  in  the  present  work  confirm,  as  has 
been  observed  experimentally,  that  the  fracture 
toughness  of  tungsten  is  affected  greatly  by  the 
type  of  impurities  segregated  to  its  grain  bound¬ 
aries.  It  should  be  noted,  however,  that  in  the 
present  work  only  cases  of  clean  grain  bound¬ 
aries  and  boundaries  fully  saturated  with  a 
single  species  are  considered.  When  two  or 
more  impurity  species  are  present  Krasko3 
showed  that  the  one  whose  energy  at  the  grain 
boundary  is  lower,  will  preferentially  segregate 
to  the  grain  boundary  and,  thus,  have  the  domi¬ 
nant  effect  on  materials  fracture  toughness.  In 
the  few  simulations  in  which  partially  contami¬ 
nated  grain  boundaries,  that  is  the  grain  bound¬ 
ary  was  not  fully  saturated  with  impurities,  were 
used,  we  observed  that  the  character  of  the 
effect  (embrittling  or  cohesion  enhancing)  ddi 
not  change  but  its  magnitude  was  lower. 
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It  should  be  noted  that  our  results  differ 
somewhat  from  the  findings  of  Krasko.3 
Specifically,  N  was  identified  by  Krasko3  as  a 
grain  boundary  embrittler  while  in  the  present 
work  this  element  was  found  to  be  a  cohesion 
enhancer.  Furthermore,  O,  A1  and  Si  were  iden¬ 
tified  by  Krasko3  as  grain  boundary  embrittlers 
while  these  elements  were  found  to  be  neutral 
in  the  present  work.  As  far  as  H  is  concerned, 
the  embedded  energy  fucntions  for  this  element 
were  not  available  and  its  effect  on  grain 
boundary  cohesion  could  not  be  studied.  We 
believe  that  the  differences  between  the  results 
obtained  by  Krasko3  and  our  results  stem  from 
the  differences  in  the  simulation  methods  used, 
Krasko3  applied  molecular  statics  and  con¬ 
strained  the  motion  of  atoms  to  the  planes 
parallel  to  the  grain  boundary.  This  constraint 
was  lifted  in  the  present  work,  and  the  simula¬ 
tions  were  carried  out  using  molecular 
dynamics.  This  method  allows  the  atoms  to 
relax  into  the  configuration  associated  with  a 
lower  energy  and  is  generally  considered  as  a 
more  appropriate  simulation  method  for  non¬ 
zero  temperatures. 


CONCLUSION 

Based  on  the  results  obtained  in  the  present 
work,  the  following  main  conclusions  can  be 
drawn: 

(1)  with  respect  to  their  effect  on  the  intrin¬ 
sic  fracture  resistance  of  the  T3  (111) 
grain  boundary  in  tungsten,  the  atomic 
species  analyzed  can  be  divided  into 
three  groups: 

(2)  Interstitially  dissolved  B,  C  and  N  are 
strong  grain  boundary  cohesion  strength¬ 
ened  which  both  enhance  intergrain 
bonding  and  give  rise  to  a  lower  level  of 
effective  loading  at  the  crack  tip  by  pro¬ 
moting  crack  tip  blunting.  There  is  some 
evidence  that  clustering  of  the  interstitial 
atoms  is,  at  least  partly,  responsible  for 
the  observed  fracture  resistance. 

(3)  O,  A1  and  Si  do  not  significantly  affect 
either  the  ideal  work  of  grain  boundary 
decohesion  or  the  crack  blunting  pro¬ 
cesses  relative  to  those  observed  in  the 
case  of  the  clean  S3  (111)  grain  bound¬ 
ary. 


(4)  P  and  S  are  strong  grain  boundary 
embrittlers  which  both  significantly  lower 
the  ideal  work  of  grain  boundary  decohe¬ 
sion  and  impede  with  crack  tip  blunting 
processes. 

(5)  Regardless  whether  the  intrinsic  fracture 
of  the  S3  (111)  grain  boundary  is  mea¬ 
sured  using  the  ideal  work  of  grain 
boundary  decohesion  or  by  the  critical 
value  of  the  F!  integral,  removal  of  dele¬ 
terious  impurities  such  as  P  or  S,  by 
either  material  purification  or  by  the  site 
competition  effect  with  B,  C  or  N  can 
have  a  great  potential  in  achieving  the 
goal  of  producing  the  ‘ductile’  tungsten. 
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